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I. INTRODUCTION 



In general, extensive systems exhibiting complex spatiotemporal dynamics including chaos may be studied as 
processes involving reaction-diffusion and convective mechanisms Q| . Analysis of the dynamics of these systems is not 
an easy task because of the large attractor dimensions involved. It would be desirable to develop ways of studying 
spatiotemporal systems using reduced model descriptions in conjunction with subsystem dynamics especially when it is 
known that extensive scaling relationships in dynamics exist as a function of subsystem size [§-|5| . Methods developed 
for low-dimensional systems may then become applicable with the concomitant advantage of simplifying the data 
requirements for studying the spatiotemporal dynamics. In this paper, we show a profitable use of this approach for 
parameter estimation with reduced model descriptions of spatiotemporal systems and which uses subsystem data for 
the characterization. The low-dimensional models can be obtained by projecting the governing equations onto relevant 
modes obtained by Karhunen-Loeve (KL) decomposition Q along with Galerkin projection Efficient ways such as 
multiple shooting boundary value algorithms [p|-[l0| that are known to curtail error propagation for low-dimensional 
chaotic and noisy data for continuous systems [ pd| , may be then reformulated, as shown here, to estimating true 
parameter values using reduced models for the spatiotemporal dynamics. 

There has been a great deal of interest in forecasting spatiotemporal time series and model identification using 
polynomial and mixed functions |L^-|h4j| . Short-term prediction of spatiotemporal dynamics by reconstruction of local 
states |l|, and KL decomposition using empirical basis functions by training amplitude coefficients using genetic 
algorithms for optimization pfj[ | have been studied and assessed. Rather than phase space reconstruction models, 
identification using some knowledge of the system structure along with nonlinear parametric regression |f~4j| has also 
been used for modeling spatially extended systems. In this context, the attempt here is to demonstrate the potential of 
a different approach that uses a multiple shooting algorithm for parameter and state variable estimations in analyzing 
spatiotemporal behavior. 

In Sec. [il], the algorithm for spatiotemporal model parameter estimation by a combination of Karhunen-Loeve 
decomposition, Galerkin's projection and Multiple Shooting (KLGMS) methodology is presented for spatiotemporal 
systems described by coupled map lattices [[t6 17 and partial differential equations. Illustrative examples using the 
KLGMS approach are presented in Sec. Ill for a single variable coupled map lattice (CML) that possesses the basic 
reaction-diffusion and convection mechanisms that give rise to complex patterns including spatiotemporal chaos and 
convective turbulence. The use of the methodology in characterizing the model from subsystem data using limited 
number of snapshots is shown and the adaptability of the method in inhomogeneous model [ [l8| identification using 
perturbation strategies is also discussed. The formalism is then applied in Sec. IV to an autocatalytic reaction-diffusion 
system described by multivariable partial differential equations (PDE) and exhibiting spatiotemporal chaos . Here 
we develop the method further for use in stringent situations when only scalar and noisy data in a single variable 
from subsystems is available for parameter estimation. 



*This preprint published in Phys. Rev E 64 056222 (2001) 
t e-mail for correspondence: ravi@che.ncl. res. in 



1 



II. KLGMS APPROACH 



Given snapshots of spatiotemporal data, u^'(n,j), for variables i = 1, 2, • • • at discrete times (n — 1, 2, • • • , M) and 
spatial nodes (j = 1, 2, • • • , L), we may obtain the fluctuating components v^'(n,j) as 

v^(n,j) = u^(n,j)-u^(j) (1) 

where 

M 

.1/ 



M 



71 = 1 



represent temporal averages for the M number of snapshots considered. The KL decomposition assumes v^'(n,j) 
may be expanded in a separable form as 



M 

- w (n,i) = E4 i) w^ i) (i)» ( 3 ) 



k=l 



where by truncating the index k to an optimum value, say N, it becomes possible to reconstruct the v^(n,j) to a 

required accuracy. In Eq. (|3|), the a^\n) are time-dependent coefficients while the (j) are spatial basis functions, 
satisfying the orthonormality condition 

(4 l \^)=S lm (4) 

with the inner product defined as 

L 

(^'^)=E#(J)^(J)- (5) 
j'=i 

The spatial basis functions can be empirical basis functions obtained from the spatial correlation matrix Fourier 
modes, Chebyshev or other sets of orthogonal polynomials, wavelets, etc. [ffl pOpl] ]. In particular, to obtain empirical 
basis functions the KL decomposition of the data v^(n, j) [Eq. || is carried out in an optimal fashion M such that 
the obtained 4>k are eigenfunctions that maximize the associated eigenvalues = (((j>k, v) 2 )- Here (•,•) denote the 
inner product as defined in Eq. (||) while (•) represents an averaging procedure that commutes with the inner product 

0. A factor t]n = J2i=i ^i/J2k=i is a measure of the energy content for increasing mode index k and decides 
an index N < M for which the series in Eq. (|^) may be truncated. More often the spatial domain is much too 

large and estimation of eigenfunctions (p^U) from a spatial correlation matrix can become quite involved. The 
method of snapshots Q helps to get over this practical difficulty by using instead the temporal correlation matrix 
Cim = i^/M)^2f = i v (ljj) v { m jj) as the kernel and reduces the KL decomposition to solving a standard eigenvalue 
problem of the form Ci m $ rn = A$; with $ being the required set of eigenfunctions {4>^\j)}- From the <p^(j) the 
respective time-dependent coefficients ai (n) may be obtained from 

b^(n)=a^(n) = (v^(n,j),^\j)) (6) 

and a full description of the dynamics as a series expansion [Eq. (^)] becomes available. Note the introduction of 

notation b^\n) in Eq. (|^) [for a£ (n)] is to convey that the time-dependent coefficients have been obtained solely 

from transformed snapshot data uW(n,j) [Eq. (Eh]. The known values of b k (n) may, therefore, be used as time-series 
for model characterization and parameter estimation studies. 

Our next step is to obtain a simpler model of the spatiotemporal system that can be used for the purposes of 
parameter estimation. Such a reduced description may be derived by Galerkin's projection of the model in conjunction 
with the KL expansion [Eq. (|h] and we shall discuss the steps involved in this procedure separately for discrete CML 
type systems and for continuous systems modeled by PDEs. 

A CML model with discrete space index j and time n may be written in a general form as 

u« (n + 1, j) = g(u<U (n, j + j'),u^ (n, j + j'), ■ ■ ■ , M ), (7) 
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where, j' = • • • , —2, —1, 0, 1, 2, • • • and g any function with fi = {fii, /U2, ■ • • , the model pa rameters. The boundary 
conditions for Eq. (Q) are generally dependent on the example being studied and in Sec. Ill we show cases involving 



periodic and open flow (as in convection) boundary conditions. In the KL-Galerkin's method we minimize the residual 

by forcing the projection of the model on the subspace of truncated basis functions <fi k 
and obtain N < M KL-Galerkin equations for the CML model described by Eq. (Q) as 



L N 



4 l) (» + 1) = E ia{ E a * (1) (•?') + (J), • ■ ■ , m) - « W (i)] <f>f U) (8) 



i=l 



where the use of the orthonormal property of the basis functions [Eq. (|j)] allows considerable simplification. The 
initial conditions for solving the map [Eq. (||)] can be obtained by forcing the initial residual to have zero projection 

on the space of basis functions, i.e., aj^(0) = J2j=i v^(0, j)4>^ (j)- 

For continuous systems in both space x and time t, the mathematical model can be written in general form as 
PDE's with appropriate boundary conditions, viz., 

u (l) (t,x) = h{u {1) (t,x),u {2) (t,x), ■■-,[!) (9) 

and the corresponding KL-Galerkin equations by projection obtained as 



N 

4 i} (*)= / ME a l 1) (*M (1) ( a; )+" (1) $ {x)dx - (10) 
i=i 



Here again N < M equations form a simpler and reduced model incorporating system parameters \i. Note 
that the summation in Eq. (^) for discrete systems becomes an integral for continuous ones, i.e., {<f>i ,(j>m) — 
Jp 1 cj)\ l \x)(j)m (x)dx with the (•,•) now denoting the usual L 2 ([0, 1]) inner product space Q defined in spatial do- 
main < x < 1. The initial conditions for solving Eq. ( fio| ) may again be independently obtained by aj ; (0) = 

We next discuss the use of the simpler KL-Galerkin models in the time-dependent coefficients ajj. ( n ) I*- 6 -' Pi- (H) 
for CML or Eq. ( [io| ) for a continuous system] , along with the known coefficient values bjj. (n) obtained from the data 
by [Eq. (^)] for the estimation of parameters /x for the respective spatiotemporal dynamics [i.e., Eq. (]?]) or Eq. @]. We 
show this is possible by using an effective algorithm in parameter estimation for low-dimensional nonlinear dynamical 
systems, viz., the multiple shooting algorithm formulated as a multipoint boundary value problem with nonlinear 
constraints for optimization . The algorithm curtails error propagation observed in chaotic dynamics and offers 
advantages in terms of number of data points required, negating effects of noise, handling missing data situations, 
parallelization and stopping at local minima during optimization [^],[l0) . 

For the CML model [Eq. (Q)], the observations in the discrete time interval [ni, um] may be chosen to form a grid 
for M multiple shooting points at rii < n 2 < ■ ■ ■ < n\ < • • • < iim modes forming (M — 1) sets of initial value problems 
in the form of [Eq. (|J)] for each of the shooting nodes ni with 1 < / < M — 1. That is, on considering k — 1, 2, • • • , N 
spatial basis modes for N < M, we obtain for a i th variable N{M — 1) maps to be solved 

af{n+ 1) = g(a^(n),^\j),u^(j)^) af '(n«) = ^(0, i 11 ) 

for an incremental time-step n — *■ n + 1. Here, Su\l) denotes the value of the i th variable at the I th shooting-point 
for k th basis mode and initial guesses for solving N(M — 1) maps of Eq. ([ll]) are taken to be the known values 

of (I) [Eq. (||)]. It may be noted that for the system governed by PDEs and data available at shooting points 
T\ < i"2 • • • < ti < • • • < tm and that are monitored at time nAt, the corresponding set of N{M — 1) initial value 
problems may be written as 

af{t)=H{af{t),4\x),u^{x),v), af '(„) = (l). (12) 

We can construct an augmented vector of initial values s and parameters \i for either model Eq. (0) or Eq. (|9|), viz., 

z = (4^(1), 4^(2), • • • , S W(M), 4 2) (1), 4 2) (2), • • • , Ml , M2 , ... , Mp ) (13) 

and attempt to minimize a least square cost function Hi (z) of the form 
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N M 



= E EE4)^ i)(n) -^ i)(a(n),M) ] 2 ' (14) 



=1,2,— fe=l n=l 



where, is a function relating components of Q in Eq. ( |TT| ) or 7i in Eq. dl2| ) and comparing to the known b£ (n) 
with cr^ the square of the standard deviation. The minimization in Eq. ( |l4| ) is carried out subject to satisfying 

oW(n,+i)-«W(i + l)-0 (15) 

so that the trajectories in the coefficients aW(rty+i) become continuous. Alternatively stated, by identifying 

yi(z) = r(s^(l), Sq^(2), • • • , s^(M), /Ui , ^2, ■ • • ,/%>) an( i y2( z ) = o^(rij+i) — sW(Z+l) we obtain, a standard nonlinear 
minimization problem |22f of the type 

min{|| yi (z) ||i I y 2 (z)^0} (16) 

z 

where the minimization of yi(z) corresponds to minimizing the cost function Eq. (|l4| ) while that for y 2 (z) implies 
satisfying the constraints imposed by Eq. (|l5|). The minimization of Eq. (^) can be carried out by starting with 
initial guess values and iterating for z using z' 9+1 ) = z^ + Az' 9 ' where, uM) G [0, 1] are damping factors. In 
doing so corrections to the augmented vector z, wz., Az^ 9 ' are obtained by solving the linearized problem 

min { || yi (.«) + M^Az^ \\l I y 2{ ^) + A.W - 0} (17) 

z t/Z C/Z 

The above Eq. @ may be solved by a suitable nonlinear optimization technique in the optimization variables, z 

[Eq. ( |l3t for arbitrary guess values for the parameters /i and initial states s^(l) — bjj. (Z). In the coding of the above 
KLGMS approach, we have employed the successive quadratic programming algorithm ]2j| coupled with numerical 
differentiation for the sensitivity matrices. For illustration, we have retained simplicity in the cost function [Eq. JT^)] 
but more effective functionals ||2jJ may be adopted in the optimizing for the parameters /i. It is to be noted that 
the methodology also allows optimizing for the sjj, (Z) even when some values of Z>j. (Z) are initially not available and 
arising due to missing snapshot data in say some i th variable. The illustrative examples presented in Sec. and 



Sec. IV bring out further the advantages offered by KLGMS approach when snapshot data availability is limited and 



possibly noise contaminated. 



III. PARAMETER ESTIMATION USING KLGMS FOR CML 



Because of their computational simplicity, CMLs are a popular and convenient paradigm for studying fully developed 
turbulence |p^ , p7| , p5| , chaos ]2q] , and pattern formation JlJ in systems. A CML model is a discrete space-time system 
with continuous state space and studies the effects of local nonlinear reaction dynamics, the coupling arising from 
diffusion due to state space gradients as well as convective effects by asymmetric coupling [jl7). Here, we consider a 
CML involving a single spatial dimension and incorporating these mechanisms as 

u(n + l,j) = (1 - D d - D c )f(u(n,j))+DJ(u(n,j - 1)) 

+y[/(«Kj + l))+/K»J-l))], (18) 

where, u(n,j), j = 1,2, • •• ,L is the state of the variable located at site j at time n for a lattice of size L, Dd the 
nearest neighbor diffusive coupling strength, and D c denoting the asymmetric coupling constant. This being a single 
variable system we suppress the index (?) in this Section. For D c = the system represents a reaction-diffusion 
system while for D c ^ mimics one with convective effects included. We assume the reaction dynamics on the lattice 
sites is governed by the nonlinear logistic function f(u) = 1 — Fu 2 , where, F is the nonlinearity parameter. Thus, 
depending on the parameter values for F, Dd and D c , a variety of dynamical patterns may be observed in Eq. jl^ ) 
and characterized as in [jlTj . We bring out the methodology for estimating parameters for selected dynamics covering 
a broad range of complexity, viz., (a) weak chaos; (b) traveling wave; (c) fully developed chaos; and (d) convective 
turbulence. Spatiotemporal data for the different cases are obtained by evolving Eq. (|l^). All the sites are given 
random initial conditions at n — and snapshots are stored after eliminating initial transients. Cases (a,b,c) are 
evolved with periodic boundary conditions, i.e., u(n, 1) = u(n, L) while for the convective case (d) the left boundary 
is assumed fixed, i.e., u(n, 1) = 1, with the right boundary open. The gray-scale images of the spatiotemporal data 
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with the parameter values yielding the data for a lattice size of L = 60 and for M = 20 snapshots, is shown in Fig. [|. 
In studies involving subsystems, only the data corresponding to the evolution of the chosen subsystems are stored. 

We obtain a KL decomposition for the spatiotemporal data v(n, j) and Table | shows the corresponding eigenvalues 
Afc, and the energy content r/k, for the data shown in Fig. [l] (a-d). The results show that for the CML exhibiting 
weak chaos and travelling wave, a smaller number of basis modes N — 3 and N — 5, respectively, are required to 
capture and reconstruct 99% of the data. For the more complex patterns, viz., fully developed chaos and convective 
turbulence the number of basis modes significantly rise to 15 for as 99% and 19 for as 100% accuracy. 

Studies with KL-Galerkin Eq. (|J) for the CML Eq. ([l8]) using the KLGMS approach did accurately and simultane- 
ously estimate the unknown parameters (F, Dd, D c ) from a few snapshots of the data. The results of convergence for 
arbitrary and different initial guesses for the parameters shown in Table |l| for the (a) weakly chaotic, (b) travelling 
wave and (c) fully developed chaos cases. The robustness is seen when parameters were successfully estimated even 
for noisy spatiotemporal data sets [Table |l| obtained by additive noise u(n,j) — u(n,j) + r\ with Gaussian distri- 
bution noise r\ s N[0,e 2 ]. The strength of the noise level used was determined by <r n oisel '&data and chosen to be 
0.01. It may be observed that noise in the data enters through the "coefficient trajectories" bk(l) that are obtained 
by the convolutions of the fluctuating data v(n,j) — u(n,j) — u(j) with the basis functions 4> k (j) via, Eq. @. Note 



that although the functional form of the CML in the form of Eq. (|18|) is single dimensional and single variable the 
procedure may be extended to situations involving multivariable mappings Eq. (Q) and higher spatial dimensions. 
The effects of considering higher spatial dimensions do not change the methodology because the KL expansion yields 
two or three-dimensional spatial basis functions <j)(i,j, k) but the Galerkin equation still retains the mapping form of 
Eq. (^) in the time-dependent coefficients (n) . Applications of the method to systems with multivariable coupling 
and scalar data are shown in the Sec. |^ studying KLGMS for continuous time systems. For brevity the results 
obtained with CMLs on these aspects are not presented. 

The presence of scaling relationships in Lyapunov exponents as a function of subsystem size have been studied . 
For KL decomposition modes, using the spatial correlation matrix, a linear relationship in KL dimension |27|] has also 
been seen. Our studies for subsystem scaling with the temporal correlation matrix, Ci m , showed some interesting 
features. We observe that T>t = maxjiV : 7]n < /} required to capture a fraction / of the total variance showed scaling 
behavior after an optimum subsystem size before saturation. The saturation occurs either due to the dynamics being 
not complicated enough to warrant all modes to be included as a function of subsystem size or alternatively when the 
dynamics is sufficiently complex that all KL modes (limited by the number of snapshots M) are required. Therefore, 
depending on the complexity of the pattern and number of snapshots, M , an optimum subsystem size exists beyond 
that only system features can be extracted reliably. The feasibility of estimating parameters by relaxing the need 
for data from the entire spatial domain was then considered. Thus, on computing T>t for the convective CML data 
[Fig. [l](G?J] as a function of subsystem size j for L = 80 we observed that beyond j = 30 there is linear scaling and 
this determines the optimum subsystem size. For this subsystem size even with a lower number of modes(w N = 15), 
parameter values could be estimated, while for larger subsystem size all KL modes need to be considered. Figure ^(a) 
shows the subsystem data for the central 31 lattice sites and used for parameter estimation purposes for N = 15. The 
accurate convergence of the estimated paramet ers F , Dd and D c with search iterations is shown in Fig. ^(7>-d) and 
reported as the homogeneous case (a) in Table III. These studies suggest that when reliability of data is poor from 
certain regions, considerable information may be gained by using only authentic data available from other subsystems 
in the spatial domain. 

A number of real situations have inhomogeneous distribution of parameter values in space and/or slowly varying 
in time domain. Studies in this context for parameter estimation were carried out and the analysis of a simple 
example is discussed here. We evolve a CML such that the sites in the left half (i.e., 1 < j < 256) have F = 2.0 
while the right half (i.e., 257 < j < 512) evolve data with F = 1.9 for L = 512. Subsystem data from each half 
[Fig. |3j was used for parameter estimation. Since the local dynamics propagate in space, the data obtained from 
both subsystems had composite features leading to inconsistent and unreliable parameter estimates. To overcome this 
difficulty we recorded data immediately after a giving a perturbation at time n (i.e., noise of strength 0.01) to the 
variable u(n,j) and then carried out KLGMS parameter estimation for each of the subsystems (left and right). The 
results presented in Table III cases b,c show that parameter estimation is now possible. Studies were also carried out 
for situations modeling F as a slowly varying parameter in time. The need to record subsystem data at optimum time 
gaps was found necessary to monitor the slow parametric changes. In real situations, repeated parameter estimations 
at sufficient time intervals can help in establishing relationships in the nature of parametric variations and this can 
considerably aid system analysis. 
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IV. PARAMETER ESTIMATION USING KLGMS FOR REACTION-DIFFUSION SYSTEM 



A basic problem in studying spatially extended dynamical systems is the quantitative comparison of experimental 
data with models based on partial differential equations. For example, in the study of pattern forming systems, the 
theoretical models usually take the form of reaction-diffusion equations that have been studied both theoretically and 
experimentally [p] j2^j29|] . For our study of parameter estimation we shall illustrate the methodology for a prototype 
reaction-diffusion model where one chemical species grows autocatalytically on another species |j30| , [l9f . This model is 
a simplification of the model of glycolysis proposed by Selkov |n| and it follows the reaction mechanism U + 2V — > 3V; 
V — > P with a continuous supply of the reactant U and removal of product P. The model has been extensively studied 
from the point-of-view of pattern formation and comparisons with features observed in experimental data have also 
been attempted [3^] . 

The reaction-diffusion mechanism yields a two variable PDE model involving concentrations t/W(t, x), u^ 2 '(t,x) of 
U , V, respectively, and for a spatially one dimensional system, we obtain: 

^ = fl„vV 1 )(i )I )-«( 1 )(t,i)( 11 ( 2 )(t,i)f + /[i- M ( 1 )(i ) ,)] 

at 

du(2) ^ x ^ = D v V 2 u^(t,x)+u^(t,x)(uW(t,x)) 2 - [f + k]u<®(t,x). (19) 
at 

Here, D u and D v are the diffusion coefficients of species U and V, with parameters / and k related to the flow of 
reactant into the system and the kinetic rate constant. The parameters f,k form a pair of bifurcation parameters that 
may be varied to obtain a host of spatiotemporal Turing patterns for unequal diffusion coefficients of the chemical 
species as seen in Q. 

In our study, we consider the situation corresponding to system exhibiting spatiotemporal chaos Fig. ||. as studied in 
p9| . For obtaining the spatiotemporal data u^(t, x), u' 2 ' (t, x), Eq. ( |l9| ) is solved numerically with Euler discretization 
in the spatial domain, with spatial length L = 1 spanning 160 spatial sites and M = 40 snapshots are stored at a 
time step At = 0.1 and with periodic boundary conditions u^\t, 0) = u^\t, L) and u^ 2 '(t, 0) = u^ 2 '(t,L) imposed. 
The initial conditions correspond to the stationary solution u^(0,x) — 1 and u^- 2 '(0, x) = except for a few central 
sites that are given a random perturbation to break the symmetry. 

Here we will also consider situations where only scalar data in a single variable (t, x) is monitored. Because the 
data in the u^ 2 -* (t, x) is not available we need to use basis functions other than empirical. In the present study, we 
choose to exemplify KLGMS using Fourier basis functions defined as 

4>f(x) = \/2sin(27rfca;) (20) 

with temporal coefficients obtained by 

b k l) (t) = af{t) = [ L v^(t,x)^\x)dx (21) 
Jo 

and use the b^\t) as observables in evaluating the least square functional in Eq.(|lJ). For the model Eq. (|l^) the 
KL Galerkin projection equations for the time-dependent coefficients, i.e., Eq. (|l0|), for modes k = 1, 2, • • • , N can be 
written as (suppressing (x) and (t)): 

„L N N N 



^ = j [^v 2 (E4 1) 4 1) +^ ) )-(E4 1) ^ 1) +" (1) )(E4 23 4 2) +- (2) ) 2 

" fc=l k=l k=l 

N 

-/(i-EW-fi^)]^ 



k=l 



r L N N N 

Jo k=l fe=l fe=l 



(22) 



fc=i 



G 



and a reduced N < M set of ODEs solved by integrating using the initial conditions discussed for Eq. (|T(]) . 

Our studies with the set of Galerkin equations Eqs. (E2n with KLGMS for estimating system parameters using the 
spatiotemporally chaotic data (Fig. ||) showed two interesting features described below. First, accurate parameter 
estimation of the diffusion coefficients {D u , D v ) did not particularly depend on the choice of (k, f) when initial 
transient data were chosen as snapshots with the diffusion mechanism playing a significant role. It was also observed 
that similar results in (k, /) were obtained using snapshots after giving a perturbation to the system state u^(t, x) 
at any time t. The second feature was that having evaluated (D u , D v ) in the above fashion the other two parameters 
(/, k) could be successfully estimated using post-transient data. These observations suggest that diffusion rates and 
reaction rates occur at differing time-scales and clearly point to the need for suitable data sampling strategies. It may 
be noted that the values of diffusion coefficients employed here lie in typical ranges. The multiple time-scale features 
discussed above may therefore be expected to be frequently present in the dynamics of spatiotemporal systems. Any 
methodology seeking model identification would need to consider this relevant aspect for parameter estimation. 

Without any ambiguity, we discuss other features of the KLGMS with reference to evaluating / and k from monitored 
post-transient data. The KL decomposition of the data set using Fourier modes for 40 snapshots showed that a single 
Fourier basis mode could reconstruct the data snapshots accurately (> 99.8%). The results of parameter estimation 
with this single mode considered showed that accurate convergence was consistently possible even when the data was 



corrupted with noise of the order of 5% and are summarized in Table IV case a. For the present reaction-diffusion 



system we have observed that the use of basis functions with the known Fourier form allows tolerance for higher noise 
levels when compared to empirical basis functions (using correlation matrices). A more practical problem arises in 
multivariable systems when only one dynamical variable is monitored. We assume that u^ 2 '(t, x) is not monitored 
and assign initial guesses for the temporal coefficients a^(t) — 0.2 and = for the multiple shooting algorithm. 
The least square functional Eq. (14|) and equality constraints are suitably modified so as to take into account only 
terms in variables u^'(t, x). Resu ts of the study presented in Table IV case b showing accurate parameter estimation 
is again possible for both /, k although with a small decrease in noise tolerance. It may be seen that the parameter 
estimation of k present only in the equation of the PDE model Eq. (|l^) is also possible. Importantly, we have 
recovered the unmonitored variable u^(t,x) using Eq. (J3J) and estimated the values of a^(t) by multiple shooting. 

Similar to the studies using CML we attempted to evaluate parameters using subsystem data with only scalar 
variable data in u^ 2 \t,x) available. An indication of the optimum subsystem size in this study using Fourier basis 

functions was suggested on evaluating the normalized power P(n s ) = /[o^^)]^ dt as a function of the subsystem size 
n s and is shown in Fig. |. The results indicate a near saturation beyond n s = 80. The results of KLGMS carried out 
with subsystem scalar data available only in and with noise (Table ^ case c) shows that parameter estimation 
within reasonable error bounds is still possible. 



V. CONCLUSION 



The results obtained using the KLGMS show that this basic framework has the necessary robustness for parameter 
estimation for spatiotemporal dynamics. We exemplify the methodology by simultaneously estimating all parameters 
of a CML and a reaction-diffusion system. Importantly, for complex dynamics and noise in the data we show that 
accurate parameter estimates are possible even from small data samples obtained from subsystems of optimal size. 
We show ways of adapting the methodology for inhomogeneous situations when parameters vary in space and time 
and by using transient data soon after perturbing the system dynamics. The usefulness of this strategy especially 
when multiple time-scales are present in the system dynamics has been discussed. The algorithm can be extended 
to situations when only scalar data is available and has the capability to recover the dynamics of the unmonitored 
variable. The study presented here should also help in the analysis and design of experiments for spatiotemporal 
systems that are often costly and difficult to perform. 
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FIG. 1. Evolved spatiotemporal data u(n,j) for the CML (j spatial grid with L = 60; M = 20 snapshots, (a) Weak chaos 
(F = 1.73, D d = 0.4, D c = 0.0); (b) Traveling wave (F = 1.5, D d = 0.5, D c = 0.0); (c) Fully developed chaos (F = 2.0, 
D d = 0.4, D c = 0.0); (d) Convective turbulence (F = 2.0, D d = 0.4, D c = 0.3). 
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FIG. 2. Parameter estimation for convective turbulence, (a) Subsystem data for the central 31 lattice sites; (b,c,d) Simulta- 
neous convergence to parameter estimates for F, D d and D c for arbitrary initial guesses (shown as y-axis labels) as iterations 
q proceed for minimizing the least square functional. 
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FIG. 3. Data from the left (F — 2.0) and the right (F = 1.9) subsystems for the inhomogeneous CML. The vertical line at 
j = 256 marks the boundary; Other parameter values D d — 0.4, D c = 0.3. 
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FIG. 4. Spatiotemporal data for the variable (t, x) in the autocatalytic reaction-diffusion system with parameter values 
/ = 0.029, k = 0.0535, D u = 0.00002, D v = 0.0001 with spatial length L = \ spanning 160 spatial sites and M = 128 snapshots 
recorded at a time step At = 0.1 is shown. 
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5. Power P(n a ) in the first mode of the temporal coefficients, normalized to the maximum, is plotted as a function of 
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TABLE I. Significance of KL modes in CML. 



Case 




Mode no. k 


At 


Vk 


(a) 


Weakly 


1 


10.5467 


0.9408 




chaotic 


2 


0.3639 


0.9733 






3 


0.2862 


0.9988 


(b) 


Traveling 


1 


13.2207 


0.9337 




wave 


2 


0.3712 


0.9600 






5 


0.1098 


0.9957 


(c) 


Fully 


1 


3.8130 


0.2690 




chaotic 


15 


0.0661 


0.9912 






19 


0.0217 


0.9999 


(d) 


Convective 


1 


4.0071 


0.2950 




turbulence 


15 


0.0458 


0.9929 






19 


0.0106 


0.9999 



TABLE II. Parameter estimation for the CML with varying dynamics. Error bounds for arbitrary initial guesses are shown. 

Case F D d 

Jo) Weakly chaotic 1.73 ± 0.01 0.40 ± 0.02 

with noise 1.74 ± 0.03 0.39 ± 0.02 

(b) Traveling wave 1.50 ± 0.01 0.50 ± 0.01 

with noise 1.54 ± 0.03 0.53 ± 0.04 

(c) Fully chaotic 1.99 ± 0.01 0.40 ± 0.01 

with noise 2.03 ± 0.04 0.38 ± 0.03 



TABLE III. Parameter estimation from subsystem CML data for convective turbulence. Error bounds for arbitrary initial 
guesses are shown. 



Case 




F 


D d 


D c 


(a) 


Homogeneous 


1.99 ±0.02 


0.41 ±0.03 


0.32 ±0.03 


(b) 


Inhomogeneous left 


1.99 ±0.03 


0.39 ±0.02 


0.29 ±0.04 


(c) 


Inhomogeneous right 


1.88 ±0.03 


0.42 ± 0.04 


0.28 ±0.04 
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TABLE IV. Parameter estimation for the autocatalytic reaction-diffusion system. Error bounds for arbitrary initial guesses 
are shown. 



Case 


Data used 


Noise level 


/ 


k 


(a) 


u [i> (x,t) 


0.00 


0.0290 ± 0.0001 


0.0535 ± 0.0001 




?/ 2) Cr t) 


0.02 


02Q1 + 0001 


0537 + 0001 






0.05 


0.0296 ± 0.0003 


0.0540 ± 0.0003 


(b) 


u^>(x,t) 


0.00 


0.0291 ± 0.0002 


0.0538 ± 0.0003 






0.02 


0.0296 ± 0.0005 


0.0565 ± 0.0003 






0.05 


0.0303 ± 0.0008 


0.0610 ± 0.0003 


(c) 


u w {x,t) 


0.00 


0.0295 ± 0.0002 


0.0540 ± 0.0002 




0.25 < x < 0.75 


0.02 


0.0307 ± 0.0008 


0.0578 ± 0.0002 






0.05 


0.0321 ± 0.0008 


0.0614 ± 0.0003 
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